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We describe repulsively interacting Bose-Einstein condensates in spatially correlated disorder 
potentials of arbitrary dimension. The first effect of disorder is to deform the mean-field condensate. 
Secondly, the quantum excitation spectrum and condensate population are affected. By a saddle- 
point expansion of the many-body Hamiltonian around the deformed mean-field ground state, we 
derive the fundamental quadratic Hamiltonian of quantum fluctuations. Importantly, a basis is 
used such that excitations are orthogonal to the deformed condensate. Via Bogoliubov-Nambu 
perturbation theory, we compute the effective excitation dispersion, including mean free paths and 
localization lengths. Corrections to the speed of sound and average density of states are calculated, 
due to correlated disorder in arbitrary dimensions, extending to the case of weak lattice potentials. 



I. INTRODUCTION 

The intriguing interplay of Bose statistics, interaction, 
and disorder is one of the most prominent problems in 
condensed matter physics, known as the dirty boson prob- 
lem [1, 2], Experimentally, it was first studied with su- 
perfluid Helium in aerosol glasses (Vycor) [3]. Over the 
past years, several groups have loaded ultracold atoms 
into optical potentials and studied Bose-Einstein conden- 
sates (BECs) in the presence of disorder under very clean 
laboratory conditions [4-9]. 

In this paper, we study the situation where Bose statis- 
tics and interaction are the dominant effects, and the 
disorder weakly perturbs the homogeneous situation. In 
this regime, the presence of a well-populated conden- 
sate makes Bogoliubov's theory [10] the most economic 
description, because it describes quantum fluctuations 
around the best mean-field approximation of the conden- 
sate. 

What happens to a homogeneous condensate if a weak 
external potential is switched on? How are the quantum 
fluctuations affected? These questions constitute the in- 
homogeneous Bogoliubov problem, a problem of notorious 
difficulty, due to the broken translational symmetry [11]. 
Bogoliubov theories for disordered BEC have been formu- 
lated by Lee and Gunn [12], Huang and Meng [13] and 
Giorgini, Pitaevskii and Stringari [ ], complemented by 
[15-21], among others. Yet, none of the existing theories 
covers spatially correlated disorder and all dimensionali- 
ties, which come into focus after recent experimental ad- 
vances [22-26]. 

Moreover, some approaches appear questionable from 
a conceptual point of view. Indeed, the primary effect of 
an external potential is to deform the condensate itself. 
This is most obvious for cold-atom BECs in traps, where 
the condensate forms in a non-uniform spatial mode that 
results from the competition between interaction, kinetic, 
and potential energy. Therefore, it is awkward, if not 
outright inappropriate, to construct a Bogoliubov theory 
in terms of fluctuations that are still defined as deviations 
from the uniform condensate of the homogeneous case 



[13, 15, 16, 19, 21]. Instead, Bogoliubov's ansatz warrants 
to first determine the deformed condensate mode on the 
mean-field level. Accordingly, we discuss the number- 
conserving deformation of the condensate caused by a 
weak external potential in Section II A. 

The deformed condensate is the vacuum of Bogoliubov 
fluctuations, to whose description we turn in a second 
step. Great care must be taken to ensure that the excita- 
tions occur in modes that remain orthogonal to the inho- 
mogeneous ground state — even in a disordered situation 
where the ground state depends on each realization of the 
random potential. We have found it helpful to tackle this 
formidable problem by a variational saddle-point expan- 
sion, a powerful method central to the solution of many 
problems in statistical and quantum mechanics [27, 28]. 
Such an expansion yields all relevant terms in a system- 
atic manner, without the need for deciding ad hoc which 
terms are to be kept or discarded. Section II B contains 
a full account of our formulation, leading to the funda- 
mental inhomogeneous Bogoliubov Hamiltonian that is 
quadratic in the fluctuations. 

We emphasize that our approach describes both effects, 
the deformation of the ground state and the scattering of 
excitations, on the same footing and to the same order in 
the external potential. Moreover, our approach involves 
quantized fluctuations that always remain orthogonal to 
the inhomogeneous Bogoliubov vacuum. A fully analyti- 
cal description is presented up to second order in disorder 
strength. 

The excitation spectrum of any system provides pre- 
cious information about its (thermo-)dynamic properties. 
For instance, the homogeneous interacting Bose gas has 
a gapless excitation spectrum. This is consistent with 
the fact that the low-energy excitations are the Gold- 
stone modes [29] associated with the spontaneous U(l) 
symmetry breaking in the BEC phase. These low-energy 
excitations are collective in character, with many inter- 
acting particles oscillating back and forth as in a sound 
wave. And really, the dispersion relation at low energy is 
linear, its slope being the sound velocity. The sound ve- 
locity is intimately linked to numerous important quanti- 
ties like specific heat and compressibility, and moreover, 



2 



by a classical argument due to Landau, equal to the crit- 
ical velocity of superfluidity [30-32] . 

Since an external potential couples to the particle den- 
sity and does not interfere with the U(l) symmetry, it is 
not expected to induce an excitation gap. However, in- 
homogeneity should certainly affect the speed of sound, 
which is a nonlinear function of the particle density. 
Thus, it is of particular interest to predict the speed-of- 
sound correction in disordered Bose gases. But curiously, 
the state of affairs for this key quantity is far from satis- 
factory. The simplest Bogoliubov theories cannot predict 
a change in excitation dispersion at all [13, 15, 16]. More 
elaborate calculations by Giorgini et al. [14] predict a cer- 
tain positive correction for uncorrelated disorder in three 
dimensions, a result which has been exactly reproduced 
by Lopatin and Vinokur [33] and Falco et al. [19]. This is 
contradicted by Yukalov and Graham [34, 35] who report 
a decrease of the sound velocity in three dimensions, even 
in the case of uncorrelated disorder. A negative correc- 
tion is also found, in all dimensions, for spatially corre- 
lated disorder with a correlation length much longer than 
the condensate healing length [36]. 

Clearly, there is a need for a unified theory that de- 
scribes the dispersion relation of Bogoliubov excitations 
in presence of disorder with spatial correlation. In section 
III, we provide such a theory, at least perturbatively for 
weak external potentials, by applying standard diagram- 
matic Green function techniques to the inhomogeneous 
Bogoliubov Hamiltonian derived in Sec. II B. We com- 
pute the ensemble-averaged disorder correction to the 
single-excitation spectrum in general, including the elas- 
tic scattering rate, and corrections to sound velocity and 
density of states. 

In Section IV, these general results are discussed in 
greater detail, with particular emphasis on the case of 
correlated disorder. Numerous analytical results are 
found in certain limiting regions of the parameter space, 
which is spanned by condensate healing length, excitation 
wave length, and disorder correlation length. Specific re- 
sults pertaining to optical speckle potentials are collected 
in Appendix A. In passing, we recover the localization 
properties of Bogoliubov excitations in one dimension as 
described earlier by Bilas and Pavloff [17] and Lugan et 
al. [18]. We briefly connect to the case of weak lattice po- 
tentials and confirm predictions by Taylor and Zaremba 
[37] and Liang et al. [38] within our formalism. We also 
reconfirm the positive correction of the speed of sound 
by uncorrelated disorder in three dimensions [14, 19, 33]. 
It turns out, though, that this result is hardly generic, 
because in lower dimensions and for correlated disorder 
in general, one always finds a negative correction. We 
confirm our analytical predictions in one dimension by 
numerical simulations on the mean field level, as well as 
by exact numerical diagonalization of the Bogoliubov-de 
Gennes equations. 

Finally, Section V concludes and closes the paper on 
some open questions. 



II. THE INHOMOGENEOUS BOGOLIUBOV 
HAMILTONIAN 

A weakly interacting Bose gas is described by the 
(grand canonical) Hamiltonian [30-32] 
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in terms of particle annihilation and creation operators 
\J/ = \J/(r) and W = W(r), respectively, which obey the 
canonical commutator relations 



[*(r),*(r')] = [&(r),&(r')) = 0, 
[*(r),*t( r ')] =5{r-r'). 



(2) 



Atom-atom interaction is taken into account in the form 
of s-wave scattering. The s-wave scattering length a s de- 
termines the interaction parameter, g = 4irh 2 a s /m in 
three dimensions, with similar relations in quasi-two and 
quasi-one dimensional geometries. We will treat the case 
of repulsive interaction with a constant interaction pa- 
rameter g > 0. This interaction potential is a good 
approximation in the regime of low energy and dilute 
gases, where the gas parameter (na 3 ) 1 / 2 is small, i.e., 
the average particle distance n -1 / 3 is much larger than 
the scattering length a s . For dilute ultracold gases (un- 
like superfluid helium), this parameter is typically very 
small K) 1 / 2 « 0.01 [39]. 

We will work in the canonical ensemble with a fixed 
total number of particles N = J d d r(& (r)*(r)>. The 
chemical potential fj, serves as the Lagrange parameter 
that has to be adjusted accordingly, as function of ex- 
ternal control parameters. One of these external control 
fields is the inhomogeneous external potential V(r). In 
the laboratory, this typically comprises a global trapping 
potential as well as, say, optical lattices and/or disor- 
der potentials. In the following, we will concentrate on 
the situation where the global trapping potential is very 
smooth, ideally a very large box, and V(r) then describes 
the local spatial fluctuations around the homogeneous 
background. 

Below a critical temperature, the Bose gas forms a 
BEC [40] , where a macroscopically large fraction of par- 
ticles populates the ground state of the single-particle 
density matrix. In the absence of interaction, this is just 
the ground state of the potential V(r), but also in a di- 
lute interacting Bose gas a well defined condensate mode 
appears, as proven rigorously in the homogeneous case 
and three dimensions [41, 42]. Within a mean-field de- 
scription (or equivalently, Hartree-Fock theory) the con- 
densate spontaneously breaks the U(l) gauge invariance 
of the Hamiltonian (1) by settling on a global phase. In 
lower dimensions and within confining potentials, quasi- 
condensates [43] exist, whose phase coherence is not of 
truly infinite range, but can extend over large enough 
distances such that the condensate shows the tell-tale 
signatures of a phase-coherent matter wave. 
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Bogoliubov's theory [10] takes advantage of the macro- 
scopically occupied ground state of the BEC and splits 
the field operator into a mean-field condensate and quan- 
tized fluctuations: 



#(r) = $(r) + (5*(r) 



(3) 



The small parameter of this expansion is again the gas 
parameter (na^) 1 / 2 [44], and for dilute condensed atomic 
gases, Bogoliubov theory proves to be a very adequate 
description. 

Following this approach, we will first describe how 
the external potential V(r) affects the condensate mode, 
strictly within mean field. In a second step, we deter- 
mine the relevant Hamiltonian of the quantum fluctu- 
ations around this modified ground state. We empha- 
size from the outset that a consistent Bogoliubov the- 
ory requires to calculate both steps to the same order 
in V(r); otherwise one runs the risk of describing only 
half of the relevant physics. Instead of deciding ad hoc 
which terms should be kept and which not, we resort to a 
well-controlled saddle-point expansion of the many-body 
Hamiltonian around the mean-field ground state. 



A. Deformed mean-field ground state 

The mean-field approach, known as Gross-Pitaevskii 
(GP) theory [32, 45], neglects the quantum fluctuations 
and replaces the field operators by a complex field \f = 
^(r), such that the many-body Hamiltonian (1) reduces 
to the GP energy functional: 
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We wish to determine its ground state as function of 
the external potential V(r). By definition, the ground 
state ^o(r) = $(r) minimizes the energy functional (4). 
It obeys the stationarity condition 6E/6^*\ Q — 0, also 
known as the stationary GP equation 



2m 



V 2 $(r) + (.g|$(r)| 2 - /i)$(r) = -V(r)$(r). (5) 



For a stationary potential V(r), the condensate's kinetic 
energy is always minimized by choosing a fixed global 
phase, thereby ruling out superfluid flow or vortices, and 
without loss of generality we may take $(r) £ 1 in the 
following. 

In the homogeneous case V(r) — 0, the repulsive in- 
teraction spreads the density over the entire available 
volume, and n = \Q\ 2 = fi/g. In the inhomogeneous 
case, however, the condensate wave function depends via 
Eq. (5) nonlincarly on V(r). Numerically, the conden- 
sate ^(r) can be computed very efficiently, for any given 
potential V(r), by propagating the GP equation in imag- 
inary time [46]. 



What analytical tools are available? If the external 
potential and the condensate wave function vary only 
very smoothly, the first term in Eq. (5), the kinetic en- 
ergy or quantum pressure, is negligible. In this so-called 
Thomas-Fermi (TF) regime, the density profile is then 
determined by the balance of interaction and external 
potential, 



n T F{r) = -\p, - V(r)] 
9 



(6) 



where V(r) < /j, and n TF (r) = else. But in a potential 
that varies on short length scales, the kinetic energy term 
becomes relevant, and the TF result no longer suffices. 

If the external potential is small, V <C gn &s [i, the 
imprint on the condensate amplitude can be computed 
perturbatively [47]. We expand 



$(r) = $ (0) +$«(r) + <I ,(2 V) 



(7) 



around the homogeneous solution = y/n in powers 
of the small parameter F//j<1. In order to maintain a 
fixed average particle density L~ d f d d r\$(r)\ 2 = n, also 
the chemical potential is adjusted at each order, 



A 4 = A* 
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We insert these expansions into Eq. (5) and collect orders 
up to V 2 /fi 2 . Because the kinetic energy e° = h 2 k 2 /2m 
is diagonal in fc-space, solving for the "J^ 1 ** and fjP' is 
best done in momentum representation, = (fe| $) = 
L -d/2 Jd d re- ife ' r $(r) and V k = (k + k'\V\k') = 



-ik-r 



V(r). 



The first-order imprint of the potential in the conden- 
sate amplitude reads 



(i) 



_ (f - 5 k v)Vk N i/ 2 



2gn 



(9) 



As expected for linear response, the shift is directly pro- 
portional to the potential's matrix element Vk- The Kro- 
necker delta stems from the first-order shift /i' 1 ) = Vq 

that compensates the potential average, ensuring ^q 1 -* = 
as required by conservation of average particle density. 

In the denominator, the comparison between interac- 
tion gn and kinetic energy defines a characteristic length 
scale of the BEC, the healing length £ = h/y / 2mgn. Fac- 
toring out gn, one is left with 2 + e®/gn = 2 + fc 2 £ 2 in 
the denominator. This term becomes constant for long- 
range potential variations with k£ — > 0, and one recovers 
the TF imprint (6) for the density n(r) = |$(r)| 2 = 
n + n^^r), in Fourier components fi^p = -V k /g. In 
the contrary case fc£ ^> 1, this denominator suppresses 
short-scale potential variations. Indeed, the condensate 
avoids rapid variations, which cost too much kinetic en- 
ergy, and responds only to a smoothed component of the 
external potential [47]. 

Figure 1 shows a ID real-space plot of the condensate 
density deformed by a rather strong Gaussian impurity 
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Figure 1. (Color online) Condensate density n(x) deformed by 
an impurity potential V(x) = V exp(— x 2 /a 2 ) (dashed black) 
with V = 0.75gn and a = 0.8£. The numerical solution of the 
CP equation (5) [solid black, under periodic boundary condi- 
tions within the shown interval] differs significantly from the 
TF result (6) [dashed green]. Including first-order (11) and 
second-order smoothing (13) terms improves the agreement. 



To second order, one finds 



$( 2 ) = _J_ V$l 1} 
k /yl/2 Z^i h—p p 



(1 - Sko)tp + 3gn 



2gn 



,(2) 
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p p 2gn + el 



L'l 



2(l-S kQ )e a p + 8gn 
2gn + e° 



(12) 
(13) 
(14) 



In Figure 1, the results of second-order smoothing are 
practically indistinguishable from the full solution of the 
GP equation. 

We note at last that even an external potential with 
zero mean causes a negative shift of the chemical poten- 
tial, 
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(15) 



potential of width a — 0.8£. For such a small impurity, 
the full numerical solution differs greatly from the sim- 
ple TF formula (6). The first-order smoothing result (9) 
already gives much better agreement. However, we need 
to push the expansion even further, in order to obtain 
consistent second-order results later on. 

Solving for the second order imprint brings about 
terms of two different types. Viewing the GP equation 
(5) as a scattering equation for the field $ [48, 49], one 
finds first a contribution from double scattering by the 
external potential V(r) with free propagation in between. 
Secondly, there is a contribution from the interaction of 
two single-scattered amplitudes Altogether, includ- 
ing the chemical potential shift, the second-order conden- 
sate deformation reads 



$i, 2) = 
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)e° - gn 



2gn- 



(10) 



For future reference, we also write down the leading- 
order terms for the Fourier components of the condensate 
density n(r) — |$(t")| 2 , the inverse field <t(r) = n/$>(r), 
as well as the inverse density n(r) — n 2 /n(r). To linear- 



order, one has of course $^ = — &y~' as well as 



-2n 



(1 - $ko)Vk 
2gn + e° 



(11) 



using the Fourier convention rip. = L~ d j d d re~ lkr n(r) 
for n and h, in the same way as for V. Eq. (11) is the 
linear response of the condensate density to the external 
potential [14]. 



A negative chemical potential shift must occur non- 
perturbatively, as can be seen by spatially integrating 
the GP equation (5) after dividing by $(r): the positivity 
of the kinetic energy entails that the chemical potential 
shift, at fixed average particle density, must be negative 
[ ]• 

This concludes our calculation of the inhomogeneous 
GP ground state, and we turn to the fluctuations around 
this condensate. 



B. Bogoliubov Excitations 

Using the Bogoliubov ansatz (3), we expand the Hamil- 
tonian (1) in powers of S^f and 8^ . To zeroth order, we 
find the GP ground-state energy E = E[$(r)]. The lin- 
ear term vanishes, because $(r) minimizes the energy 
functional (4). The relevant contribution is then the 
quadratic part, E = Eq + H. Third-order and forth- 
order terms in the fluctuations are neglected. They de- 
scribe interaction between the excitations and become 
only relevant for larger densities or higher temperatures 

[ ]■ 

For reasons that will become clear in Sec. II B 5 
below, the inhomogeneous Bogoliubov Hamiltonian is 
best expressed in density-phase variables. From ^/ = 
exp{i5(p}y/n + Sn = $ + 5h/2Q + iQStp + . . . follows 

5h(r) = $(r){s&(r) + 6$(r)} , (16a) 

Mr)= 2$^){^ t(r)_ ^ (r) }' (16b) 

up to higher orders in 5^/. The commutators (2) imply 
that density and phase (fluctuation) operators are con- 
jugate, [Sn(r),5<p(r')\ = i8(r — r'). By expanding the 
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many-body Hamiltonian to second order in the fluctua- 
tions around the mean-field solution, the relevant Hamil- 
tonian for the excitations is found as [51] 



H = 



, 5n \ 2 [V 2 $(r)] 
2m \ ' 2$(r) J + 4$ 3 (r) 

$ 2 (r)(V<5<^) 2 



Sn 2 



9 r~2 

- on 



(17) 



With Eq. (17), the problem is reduced to a Hamiltonian 
that is quadratic in the excitations. To this order, there 
are no mixed terms of Sn and 6<p. The perturbing po- 
tential V(r) does not appear directly. Instead, it en- 
ters nonlinearly via the condensate function 3>(r), which 
can be predetermined by solving the GP equation (5) or 
calculated perturbatively, as explained in the previous 
Sec. II A. 

Before further discussing the impact of the external 
potential, we briefly consider the excitations of the ho- 
mogeneous system. 
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Figure 2. Bogoliubov excitations in the homogeneous conden- 
sate, (a) Bogoliubov dispersion relation (24) (solid), (b) Den- 
sity of states (26) in three dimensions (solid). In both panels, 
dashed and dotted lines show the low-energy and high-energy 
asymptotics, respectively. The gray dashed lines foreshadow 
the disorder corrections provided in Sec. IV. 



1. Homogeneous Bogoliubov Hamiltonian 

In the homogeneous case V(r) — 0, the Bo- 
goliubov Hamiltonian (17) becomes translation invari- 
ant and thus diagonal in the momentum represen- 
tation Shk = L~ d l 2 J d d re~ lk ' r 5h(r) and S(p k = 
L -d/2 jd d re- lk r Stp(r): 



£(0) = ^ 



n r> t r , 2gn + e? + 
ne k Sip k Sp k + — — — Sn k Sn k 



(18) 



with e k = h 2 k 2 /2m. This Hamiltonian looks diagonal, 
but the Heisenberg equations of motion for Sn k = Sn_ k 
and 5cp k = <^-fcj which obey 



[5h k ,§<p k ,] = iS, 



kk' i 



(19) 



are still coupled. This is resolved by a Bogoliubov trans- 
formation [10], coupling density and phase fluctuations 
to quasiparticle creation and annihilation operators 7^ 
and %: 



Ik 

7i, 



= A k 



i\fn SCp k 
Sh k /2y/n 



A k = ak % 
-a k a k 



(20) 

A transformation of this kind, with the free parameter 
a k , guarantees that the quasiparticles obey bosonic com- 
mutation relations 



;,t 



7fe'7fe' 



|7fc ; 7fc'J = S k k>, [7fc:7fe'] 
The Hamiltonian (18) becomes diagonal 



# (0) = E ^7*, 



0. (21) 



(22) 



by choosing 



a k = 



1/2 



2 + k 2 ^ 



(23) 



The excitations are found to have the famous Bogoliubov 
dispersion relation [10] 



e k {2gn + e k ) = gnkZs/2 + k 2 e, (24) 



plotted in Figure 2(a) for reference. 

In the high-energy or large- momentum regime kt; 3> 1, 
the excitations are essentially free particles with disper- 
sion e°, shifted by the condensate background energy, 



efe « £fe + gn. 



(25) 



In the low-energy regime, the interaction dominates over 
the bare kinetic energy. A single excitation involves many 
individual particles, comparable to a classical sound 
wave. Indeed, the dispersion relation is linear, e k — hck, 
with the bare sound velocity c = \fg~nfm. According 
to a classical argument due to Landau [31, chapter 10.1], 
this linear dispersion at low energies implies superfluidity 
with critical velocity v c — mm k (e k /hk) — c. 

The transition from sound-wave to single-particle ex- 
citations also reflects in the density of states (DOS) 



Sdkt 1 
{2ir) d 



Ok 



6(6) 



(26) 



with Sa — 2, 2-7T, 47r the surface of the d-dimensional 
sphere in d = 1,2,3, respectively. The k- vector at energy 
e is given by k 2 ^ 2 = [1 + (e/gn) 2 ] 1 ^ 2 — 1. Eq. (26) shows 
the transition from a sound-wave DOS p S w(e) oc e 



d-l 



particle DOS p p t(e + gn) oc £2 x 5 
ure 2(b) for d = 3. 



to 



as illustrated in Fig- 
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2. Inhomogeneous Bogoliubov Hamiltonian 

Let us come back to the Hamiltonian (17) including 
the full imprint of V(r) in the condensate $(r). This in- 
homogeneity breaks translation invariance, so the Hamil- 
tonian cannot be diagonal in momentum representation. 
However, it is still quadratic in the fluctuations without 
any term mixing Sn and Sip, with the general structure 



H = 



2 XI \ nS< Pk S kk " S r k ' 



k,k' 



1 

4n 



5rV k R kk ,5h y \. (27) 



The coupling matrices Skk' and R kk > contain the rele- 
vant information about the Fourier components of the 
condensate field $(r) and its inverse <t(r) — n/Q(r), as 
well as of their gradients. 

There is only a single term involving the phase gradi- 
ents in the Hamiltonian (17), proportional to the density. 
Upon Fourier transformation, the coupling matrix reads 



Skk' — k • k Tlk—k' 

run 



(28) 
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at 



Figure 3. Universal Bogoliubov scattering vertex (31). Bo- 
goliubov excitations (7fc,7l fc ) are scattered by an effective 
vertex, nonperturbatively determined by the deformed con- 
densate and its inverse, ($ q ,$q). 
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Y kk > 



Y kk > 

w kk , 



fcfe' 



(31) 

depicted in Fig. 3. This brings the inhomogeneous Bo- 
goliubov Hamiltonian (17) in the form H = + 



Its diagonal elements are S kk = 2e°, to all orders of V(r), 
by conservation of average density. In the homogeneous 



case, it reduces to S 



(0) 

kk' 



2e®5 kk <. As a function of the 



condensate field components, it can be rewritten as 



Skk' — 



2ff v- 
L d 2^ 



k ■ k' £ 2 <Z> k -p$p-k' 



(29) 



In the Hamiltonian (17), the density fluctuation 8n 
appears in several, complicated looking terms. But the 
corresponding coupling matrix R kk > can be brought in a 
form very similar to (29) by using the Fourier components 
of the inverse field: 



Rkk' ^ ^ T'kpk' ^k—p^p—k' 

P 



AgnS, 



kk' 



(30) 



With f kp y = [p 2 + 2(fc' - p) ■ (k - p) + \{k' - p) 2 + 

|(fe — p) 2 ]£ 2 - In the homogeneous case, it reduces to 

Ruli = 2(2gn + t\)8 kk < . In the inhomogeneous case, the 
background-mediated coupling between fluctuations, as 
expressed by Eqs. (28) and (30), is non-perturbative in 
the potential strength [as long as the Bogoliubov ansatz 
(3) is valid]. These expressions hold for arbitrary po- 
tentials, if only the condensate $ and its inverse $ are 
correctly determined. Notably, the commutation rela- 
tions (19) for the fluctuation operators, as defined in Eq. 
(16), remain valid in the inhomogeneous setting. 

For further analysis in terms of Bogoliubov quasiparti- 
cles, we transform to the Bogoliubov basis (20) of the ho- 
mogeneous case. We separate the homogeneous contribu- 
tion from the inhomogeneous contribution in the coupling 



matrices, S kk ' = S 



(0) 
kk' 



Sjfij and R kk > 



R 



(0) 
kk' 



R 



(V) 
kk' 



and define the effective Bogoliubov excitation scattering 



E 

k 



H = > ekl k lk- 



2 <H 

k.k' 



(7fe>7-fe) 



W 
Y 



W 



Ik' 

kk' V i-k' 



(32) 

Let us reflect on what has been achieved at this point. 
By a saddle-point expansion of the general many-body 
Hamiltonian (1), we have derived the Hamiltonian de- 
scribing the dynamics of Bogoliubov excitations in inho- 
mogeneous external potentials. These excitations are de- 
fined as in the homogeneous case, with now an inhomoge- 
neous contribution to the Hamiltonian, the coupling ma- 
trix Vfcfe', that provokes scattering between different k- 
modes. This coupling has a much richer structure than a 
simple potential scattering term V q a k+q a k for single par- 
ticles. It is both nonlinear in the potential and contains 
off-diagonal contributions, because the underlying con- 
densate background depends nonlinearly on the potential 
and mediates anomalous scattering between quasiparticle 
excitations. Nevertheless, we have managed to identify 
the relevant scattering vertex V, which allows us to set up 
a systematic perturbation theory. As a first step for fully 
analytical calculations, we have to expand the scattering 
vertex to lowest orders in V. 



3. Perturbative expansion of the Bogoliubov scattering 
vertex 

We expand the scattering matrix elements in powers 
of the inhomogeneous potential V kl using the smoothing 
theory exposed in Sec. II A: 



* * 

V= ® 



V 1 



•V 



(2) 



(33) 
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Each dashed dangling line here represents the bare ex- 
ternal potential V q [not to be confounded with the dash- 
dotted double lines representing the background conden- 
sate fields ($q,<i>q) in Fig. 3]. The scattering vertices 
of order j can be derived systematically by 

(i) computing the ground state (7) to the desired or- 
der; 

(ii) computing the inverse field $ q = (n/(fr) q ; 

(hi) collecting all terms of order j in Eqs. (28) and (30) 
in order to obtain S^' and RP', and 

(iv) applying the transformation (31) in order to obtain 
and 

Let us make this procedure explicit for the first orders 

i = i,2. 

The first-order scattering amplitudes 



kk' v k-k' 



Y$=y ( klv k - k ', (34) 



are both directly proportional to the potential's matrix 
element V k _ k > , as required by conservation of momen- 
tum. All information about the interaction and the back- 
ground condensate is factorized into the envelope func- 
tions 



( (i) _ (1 - S kk ')a k ak>i 2 
,J kk> - 2 + e(k' -k) 2 



k 2 



fc' 2 kk' - 



k ■ k' 



2„2 



(35a) 



(l) _ (1 - S kk ')a k ak'^ 2 
2 + £, 2 (k' -fc) 2 



y kk 



fc' 2 kk' 



k ■ k' 



2^2 



a, .a 



k u k> 



(35b) 



with afc from Eq. (23). We have previously studied the 
scattering of Bogoliubov excitations by an isolated im- 
purity, as described by these matrix elements [E ]. An 
interesting feature of the transition from sound-like to 
particle-like excitations is that the amplitude for elastic 
scattering |fe| = |fc'| by an angle 9 is proportional to the 
remarkably simple envelope function 



w 



(i) 

kk' 



k'=k 



4 k 2 e(i 



cos6> 



fc 2 £ 2 (l 



(1-<W)- (36) 



=: A{k£, 



This envelope is responsible for a node in the scattering 
amplitude at cos6>o = fc 2 £ 2 /(l + fc 2 £ 2 ), thus interpolating 
between the p-wave scattering of a sound wave with 9q = 
7r/2 and the s-wave scattering of a single particle [51]. 
The second-order couplings 



^fefe' — „„ ^ S kvk'Vk-pV p -k' 
9 n p 



R 



(2) 
kk' 



1 X! ' » r " 



gn 



kpk' V k-p Vp-k' 



(37) 
(38) 



feature the kernels 



M 

kpk 

r (2) 
kpk' 



, = 2£ 2 fc • fc' 



{(k-k') 2 + (k~p) 2 + ( P -k') 2 }e 

[2 + (fc - k') 2 e] [2 + (fc - p) 2 e] [2 + (p - fc') 2 ^ 2 ] 

2e 2 |p 2 + 2(fc - P ) ■ (fc' - P ) + - pf 



(1 - <W)(! _ S kp )(l - 5 pk <), 



+ 2(fc 2 + fc' 2 -fc-fc') 



3 + (fc - k'fe - e(k - P ) 2 (i - s kk >) \ i-s, 



2 + (fc-fc') 2 £ 2 



^kp 



i-s, 



pk' 



2 + (k- P ) 2 e 2+{p-k') 2 ^ 2 ' 



(39) 



(40) 



Later, the ensemble average over the disorder will restore 
translation invariance. The relevant diagonal elements 



are s 



(2) _ 
kpk 



(as required by Eq. (28) and conservation 



of particle number) and 



(2) 
kpk 



2£ 



2 p 2 + 3(k-p) 



3fc z 



[2+(k- P ) 2 e] 2 



(1 - S kp ). 



(41) 



Finally, the matrix elements are transformed according 
to Eq. (31), which yields the second-order diagonal scat- 
tering amplitudes 



W, 



(2) 



kk 



Y, 



(2) 



kk 



E 



w kik v k- 



pVp- 



(42) 



where u4pfe — a 2 r k ^ k / 'Agn in terms of Eqs. (23) and (41). 



,2j2) 



4- One- dimensional setting 

Before proceeding with the general theory, we briefly 
digress into one dimension to discuss the link to previous 
works. Using the general formalism developed so far, we 
can easily compute the reflection of a Bogoliubov excita- 
tion by a delta-like impurity V(x) = Vuq5{x) to lowest 
order in Vctq. In one dimension only exact backward 
scattering occurs, and the problem involves the coupling 
element W— k ,k- In the Born approximation, to second 



order in V, we find the transmission 



T = 1 — 



r 2 



4g 2 n 2 (/fc 2 e 2 + l) 2 



(43) 



The impurity becomes perfectly transparent at long wave 
lengths fcoo — > 0. At the crossover k£ sw 1 from sound 
waves to particles, strong backscattering leads to a trans- 
mission minimum, whereas at high energies the trans- 
mission increases again. The result (43) is in quantita- 
tive agreement with the results of Bilas and Pavloff [17], 
taken in the limit of a weak impurity. It may be useful 
to note that within our formalism, we do not require the 
explicit knowledge of Bogoliubov eigenstates, nor distin- 
guish between propagating and evanescent modes; all the 
physics is built into the effective scattering vertex V. 

Also Kagan, Kovrizhin, and Maksimov [ ] have con- 
sidered tunneling across an impurity, which in their case 
suppressed the condensate density very strongly. We 
agree in the aspect of perfect transmission at low en- 
ergies. At high energies, however, Kagan et al. do not 
find a revival of transmission. This is reasonable because 
their strong impurity deeply depresses the condensate on 
the spatial scale of £, and for wave lengths shorter than 
£, transmission remains suppressed. 



5. Appropriate basis for the inhomogeneous Bogoliubov 
problem 

Before deriving physical quantities from the effective 
Hamiltonian (32), which will be the subject of the follow- 
ing section III, it remains to justify our choice of basis 
for inhomogeneous Bogoliubov excitations. 

The Hamiltonian (17) is quadratic in the fluctua- 
tions, no matter whether written in terms of the single- 
particle basis S^>(r), S^f'(r) or the hydrodynamic ba- 
sis Sh(r),5<p(r), and thus can always be diagonalized: 
H = ^ u huj v j3ll3 v . Here, the eigenmodes v are popu- 
lated by bosonic quasiparticles, for which 



frv = J d d r 

[flfj,: ft/] — fifivi 



u*(r)<5*(r) -I- w*(r)<5* t (r) 
01 Jt] = [4,/3„]=0. 



(44) 
(45) 



The eigenfunctions u u (r) and v u (r) are solutions of the 
Bogoliubov-de Gennes equation, a non-Hcrmitian eigen- 
value problem: 



H(r) gn(r) 
gn(r) H(r) 



0-3 



thbJ,, 



u v {r) 
v„(r) 



0, 



(46) 



|^ V 2 + V(r) -fi + 2gn(r) and the Pauli 
In the case of broken translation 



with H{r) = 
matrix 03 = ( — 1 ) 
symmetry, these modes are not indexed by a wave vector 
fc. They do, however, fulfill the bi-orthogonality relation 

/ d d r K(r) Ufl (r) - i£(f>„(r)] = 5^ (47) 



and the orthogonality with respect to the condensate [53] 
d d r$*(r) [u v (r) - v v {r)} = 0. (48) 



This latter relation expresses the bi-orthogonality (47) 
with respect to the zero-frequency Goldstone mode re- 
lated to the spontaneously broken C/(l) symmetry of the 
BEC [29, 54]. The orthogonality relations allow the in- 
version of Eq. (44) : 



(49) 



In presence of external inhomogeneities, and in partic- 
ular for a disorder potential to which we will turn shortly 
[Sec. Ill], this eigenbasis explicitly depends on each po- 
tential realization, which renders it useless for analytical 
calculations. Instead, we construct a basis starting from 
the plane waves that diagonalize the clean Hamiltonian, 
while satisfying the orthogonality relations (47) and (48), 
even in the inhomogeneous case. The price to pay for 
using a plane-wave index is of course that the disorder 
leads to scattering between these eigenstates. But this 
is always the case in disordered systems, and standard 
perturbation theory applies. 

Still, one has essentially two choices, (i) One can de- 
fine Bogoliubov operators by an expansion over single- 
particle plane-wave modes: 



^w = EK ) (^ ) -4 0) w^ )t ) 



with 

„.(°)/ 



Uk 



L -d/2 e ik-r 



(0) 



(r) = v k 



j —d/2 ik-r 

•J-J fc- 



(50) 



(51) 



The coefficients = ^(a^ 1 + a k ) and v k — ^(a^ 1 — a/-) 
are custom-tailored to satisfy the bi-orthogonality (47) 
for all k 0, because n\—v\ = 1. But if now the disorder 
potential is switched on, the condensate $(r) is deformed 
and not orthogonal to the plane waves anymore. Testing 
the condition (48), we find 



A d r^{r)\uf(v)~vf{v)\ = (u k -v k )$_ k jt 0. (52) 



This overlap with the ground state has disastrous conse- 
quences for the theory. If one tries to work with these 
operators, the coupling W^k 1 diverges for k — > 0, and per- 
turbation theory will break down, no matter how small 
the external potential. 

(ii) One can define the excitations via the hydrody- 
namic fluctuations (16), 



6V(r) = 



Sn(r) 
2$(r) 



i$(r)5ip(r). 



(53) 



Contrary to case (i), now the disorder is present from 
the outset, such that the fluctuations originate from the 
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disorder-shifted reference point &(r), which corresponds 
to the Bogoliubov vacuum of the true excitations (44). 
The inverse Bogoliubov transformation (20), 



Green functions 



5h k = a k y/n (% + 7-fe) > 



2iak\/n V K J 
then leads to a decomposition of the form (49), 



over mode functions 

I ( $(r) a k ^\ e ik r 

Uk(V) ~ 2 \a k ^i W) J L d/2 '' 

, . 1 / $(r) ak\/n\ e tk r 
v k {r) = ' 1 



2 Va/cV" / L d / 2 



(54) 
(55) 

(56) 

(57a) 
(57b) 



In the homogeneous case $ = i/n, these functions re- 
duce exactly to the plane- wave amplitudes (51). In the 
inhomogeneous case, the plane waves are found to be 
modified in such a way that the modes still satisfy the 
bi-orthogonality (47). Moreover, they also respect the or- 
thogonality to the deformed ground state (48), because 
$(t) [u k (r) — v k (r)] is a plane wave with zero spatial av- 
erage for all k 0. 

In conclusion, the Bogoliubov quasiparticles defined in 
terms of density and phase via (54) and (55), or equiva- 
lently by (20), fulfill all requirements for the study of the 
disordered Bogoliubov problem. They can be labeled by 
a wave vector k, which is independent of the disorder re- 
alization V(r), they fulfill the required bi-orthogonality 
relation (47), and most importantly, they decouple from 
the inhomogeneous condensate ground state. 



III. MODIFIED EXCITATION DISPERSION 

In the previous section, we have set up the general for- 
malism for describing Bogoliubov excitations in a weak 
external potential, by deriving the relevant Hamiltonian 
(32) in the form H = + H (v \ where ffW describes 
the clean system, and the disorder. This structure 

permits using the machinery of perturbation theory [55- 
57]. Presently, we explore the consequences for the exci- 
tation dispersion relation and the corresponding density 
of states. These quantities can be computed via zero- 
temperature single-excitation Green functions, for which 
we calculate the self-energy to order V 2 . From the self- 
energy, we determine physical quantities like mean free 
paths and corrections to the speed of sound. Since we 
will mainly focus on the case where V(r) is a disorder 
potential, we use a notation adapted to that scenario 
in the following. But the general theory applies to ar- 
bitrary potentials and notably covers the case of weak 
lattice potentials, to which we devote a brief discussion 
in Sec. IV B 4 below. 



The matrix structure of the scattering vertex V defined 
in Eq. (31) suggests introducing the Bogoliubov-Nambu 
(BN) pseudo spinors t k = (7fe,7_ fe ) t in terms of which 
the Hamiltonian (32) takes a more compact appearance: 



H = 



fefe' 1 k' 



(58) 



k,k' 



The Heisenberg equation of motion for the BN spinor 
reads 

d - 

ih gl^k = 0-3 £ [efc<W + Vfefe'] t k ' ■ (59) 



k' 



The multiplication by the Pauli matrix 03 is character- 
istic for the dynamics within the Bogoliubov-de Gcnncs 
symmetry class that describes bosonic excitations of in- 
teracting systems [58]. We see that the Bogoliubov ex- 
citation in mode k is scattered to mode k' by the po- 
tential Vfcfc', the momentum transfer being provided by 
the underlying condensate <I> and its inverse profile $, 
as represented by the vertex in Figure 3. The effective 
excitation spectrum belonging to the equation of motion 
(59) can be derived by studying the corresponding Green 
function. 

Many-body Green functions contain all the informa- 
tion about how a quasiparticle created in state k' at time 
propagates to state k where it is destroyed at time t. 
For the present, we need only the retarded Green func- 
tions at temperature T — 0. Taking advantage of the 
Nambu structure, one defines a matrix-valued Nambu- 
Green function [55] 

«-w-^w,fS,(P)]>-(^^w 

(60) 

from the single-(quasi)particle retarded Green function 



1 



G k y{t) = -Q{t)( [%{t),%,{$)]) (61) 



and the anomalous Green function 



1 



F kh >(t) = Jh e(t)([^ k (t),f h ,(0)]). (62) 

Here, (•) stands for the expectation value in the Bogoli- 
ubov vacuum |0) defined by ^ k |0) = for all k. The 
equation of motion of Q under the Hamiltonian (58) reads 



ih-g = a z 5{t) + g 3 [e + V]G. 



(63) 



In this compact notation, e kk > = £ k 5 kk >t. 

B. Perturbation theory 

In absence of disorder V = 0, the equation of mo- 
tion (63) is readily solved in frequency domain. The 
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anomalous Green function F^,(lu) — vanishes, and 
the conventional retarded Green function is found as 

t(0) 



G L'M = S kk>Gok(u) with 



GofcM = lim 



1 



1 



?)->o — £fe + zt; fej — efe + iO 



The infinitesimal shift -MO stems from the causality fac- 
tor 0(f) that is characteristic for the retarded Green func- 
tion (61). The Nambu-Green function for the clean sys- 
tem thus reads 



The general structure of the first contribution VWCJqVW 
is well-known from single particles in disorder [59—62]: 
the particle is scattered once by the bare disorder into a 
different mode and then scattered once more, back into 
(64) the original mode. The second contribution V( 2 ) is spe- 
cific to the Bogoliubov problem and the nonlinear back- 
ground of the GP equation (5): it describes the single 
scattering of the excitation by a background fluctuation 
that is itself second order in the disorder amplitude. 



( G ak (oj) 



(65) 



V Go*(-w) 
With this, equation (63) can be written in the form 

g=[Gv 1 -Vy 1 (66) 

which is a suitable starting point for diagrammatic per- 
turbation theory. 

Equation (66) permits a series expansion in powers of 
V for the full Green function, ensemble-averaged over the 
disorder: 



G = Go + GoVGo + GoVGoVGo 



(67) 



Without loss of generality, we assume in the following 
that the disorder potential is centered, i.e. V = 0. Then, 
second-order and higher moments of the disorder poten- 
tial have to be computed, V^V^, Vfc 1 Vfc 2 Vfe 3 , etc. De- 
pending on the disorder distribution, these moments may 
factorize into independent terms. E.g., the moments of a 
Gaussian random process factorize completely into prod- 
ucts of pair correlations. Thus, the series (67) contains 
reducible contributions from products of disorder corre- 
lations that can be separated into independent factors 
by removing a single Green function Go- This redun- 
dancy can be avoided by defining the self-energy £ via 
the Dyson equation 



G = Go + Go^G- 



(68) 



The self-energy contains precisely all irreducible con- 
tributions of the disorder-averaged right hand side of 
Eq. (67). Moreover, it directly describes the disorder- 
induced corrections to the spectrum, as becomes evident 
from the formal solution G = Gq 1 — S. 

In principle, any desired order in the disorder poten- 
tial V of the series £ = £^ + E^ 2 )_+ E< 3 > + . . . can be 
determined by first expanding £ = V + VGoV + ■ • ■ into 
powers of the nonlinear scattering vertex V and then us- 
ing the perturbative expansion (33), while retaining only 
the irreducible contributions. In practice, of course, the 
number of diagrams grows very rapidly with the order. 
Since all first-order terms vanish by virtue of V = 0, we 
concentrate on terms of order V 2 . This so-called Born 
truncation of the full series is valid for weak potentials. 
We find two contributions: 



C. Self-energy 

(2) 

We focus now on the upper left-hand block £]y of the 
Nambu self-energy matrix, which relates to the normal 
retarded Green function and describes the change in the 
quasiparticle dispersion. 

Spelling out the two contributions (69) in terms of the 
first- and second-order scattering matrix elements (34) 
and (42), we find 



,(2) 
J llfefc 



(w)v r fc_ p y p _ fc '. 



(70) 



Under the sum, we distinguish two essential factors, the 
potential correlator and a kernel function. Concerning 
the potential, we assume for notational convenience that 
the disorder is homogeneous and isotropic under the en- 
semble average, with a fc-space pair correlator 



V q V- q > = L- d 5 qq ,V z a d C d (qcj) 



2_d/ 



(71) 



The dimcnsionless function C'd{qcr) characterizes the po- 
tential correlations persisting on the length er. Our 
formulation allows for a straightforward extension to 
anisotropic disorder [26] or lattice potentials, see Sec. 
IV B 4 below. 

Because this disorder average restores homogeneity, we 
only need the kernel function for k' = k: 



Zkpk (w) 



(1)12 
kp J 



+W W 

huj-e v + i0 hu + e v + i0 kpk ' 



(72) 



with the envelope functions defined in Eqs. (35a), (35b), 
and (42). This kernel depends solely on the heal- 
ing length £. Finally, the retarded normal self-energy 



,(2) 
J llfcfc 



,(lu) = (j) takes the functional form 



£(jfe,u;) = V 2 a d 



d d q 
(2tt 



d Z k (k+ q )k{u)C d {q(j). (73) 



This expression is the second main achievement of the 
present work. All physical results presented below follow 
by straightforward calculations from here. 



D. Disorder-modified dispersion relation 



£( 2 ) = VM^oVW + V( 2 ) 



In order to grasp the significance of the self-energy, 
(69) it is useful to define the spectral function S(k, oS) = 
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— 2lmG(k,uj) [ ], which contains all information about 
the frequency and lifetime of the excitations. In the clean 
system with dispersion e k = hcu, the spectral function is 



So(k, lo) = — 2ImGo(fc, lo) = 2ir5(huj — e k ). 



(74) 



In presence of disorder, this function gets modified, while 
retaining its normalization J(dhLu/2ir)S(k,uj) — 1; this 
allows interpreting the spectral function as the energy 
distribution of a quasiparticle with wave vector k. To 
leading order in V 2 /(gn) 2 , 



S(k,co) 



-2ImS(fc,w) 



[huj - (e k +Re£(k,uj))Y + [Im£(fc,u)]' 



(75) 



The self-energy's real and imaginary part enter in charac- 
teristic ways. First, the Bogoliubov modes are expressed, 
according to the reasoning of Sec. II B 5, in the plane- 
wave basis, which is not the eigenbasis of the Bogoliubov 
Hamiltonian in presence of disorder. Thus, k are not 
"good quantum numbers" , and the Bogoliubov modes 
suffer scattering. This broadens their dispersion relation, 
or equivalently implies the existence of an elastic scat- 
tering rate (inverse lifetime) 7fe = r^T 1 = — 2Im£(fc)/fi.. 
Here, the notation £(fc) = £(fe, e k ) indicates that one 
can take the self-energy on shell to the lowest order V 2 
considered. In terms of length scales, the scattering rate 
defines the elastic scattering mean free path l s via 



17 1 = — 



lk _ -2ImS(fc) 

<>!, ' I, 



(76) 



with the usual definition of the group velocity, Hv g = 

d k e k = 2. 9 <(i + k 2 e)[2 + k 2 er i/2 . 

Secondly, the peak of the spectral density is shifted to 



e k = e k +ReE(fc), 



(77) 



which defines the disorder-modified dispersion relation, 
again to order V 2 . Notably, Eq. (77) describes the im- 
pact of disorder on the speed of sound in the low-energy 
regime. Indeed, a short calculation shows that the ker- 
nel function (72) behaves like Zk^+ q )k ~ k as k — > 0, 
such that there is always a finite correction to the speed 
of sound. Within our theory, the disordered potential 
conserves the linear character of the dispersion relation 
at low energy, as required by the existence of the zero- 
frequency Goldstone mode due to the spontaneously bro- 
ken U(l) symmetry. 



E. Average density of states 

The quasiparticle dispersion enters practically all ther- 
modynamic quantities that determine how the disordered 
condensate responds to external excitations, both at zero 
and finite temperature. Often, one only needs to know 
the density of states. In a disordered system, the spec- 
tral function (75) — remember its role as the probability 



density for a Bogoliubov quasiparticle k to have energy 
Hu> — determines the average density of states (AVDOS) 
per unit volume as [55, 57] 



p(hu>) = 



d d k S(k,ui) 
(2n) d 2tt 



(78) 



As function of frequency and for weak disorder, the 
spectral function (75) is very well approximated by a 
Lorentzian centered at e k /h with small width j k <C e k - 
The following Sec. IV will show that the relative scat- 
tering rate jk/ek of low-energy, sound-wave excitations 
tends to zero, and that the main effect of disorder is 
the dispersion shift (77) (in contrast to the case of sin- 
gle particles in disorder, where the scattering rate is the 
dominant quantity [62]). In the sound-wave regime, we 
can therefore approximate S(k,u>) = 2i:d(huj — e k ) in Eq. 
(78). With this, the shift p(e) = p(t) + Ap(e) from the 
clean DOS in d dimensions [Eq. (26)] reads 



Ap(e) 



d+k lk_ 


ReS(fc) 




kd k e k 



(79) 



k=k e 



and is thus found to be a function of the dispersion shift 
(77). 



F. Parameter space of the disordered Bogoliubov 
problem 

Because the self-energy is evaluated to order V 2 , all 
the corrections it implies will be of this same order. For 
notational brevity, we henceforth denote the small para- 
meter of this expansion by 



r 2 



(gn)'- 



< l. 



(80) 



Furthermore, we observe that the self-energy (73) de- 
pends on three different length scales: the excitation 
wave length A = 2ir/k, the healing length £, and the 
disorder correlation length a. The resulting physics can 
only depend on the value of these lengths relative to each 
other: 

• The correlation parameter £ = cr/£ indicates 
whether the disordered condensate background is 
in the Thomas- Fermi regime (( > 1) or in the 
smoothing regime (( -C 1), see Sec. II A. 

• The reduced wavenumber k£ indicates whether the 
excitations are sound waves (fc£ <C 1) or particles 
(fcf» 1), see Sec. IIB1. 

• The parameter ka discriminates the effectively 5- 
correlated regime (ka <C 1) from a very smooth 
scattering potential (ka ^> 1). 

These parameters are not independent, for any given 
two of them determine the third, e.g. (ka)/(k^) = <r/£. 
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Figure 4. (Color online) Parameter space of the disordered 
Bogoliubov problem, spanned by three length scales: healing 
length f, excitation wave vector k, and disorder correlation 
length a. At opposing vertices, the three dimensionless para- 
meters fc£, ka, and <r/£ take their extreme values or oo. On 
the six edges, one of the length scales itself is either or oo. 




0.5 1 1.5 2 2.5 3 3.5 4 

Figure 5. (Color online) Inverse elastic scattering mean free 
path (81) as function of Bogoliubov wave number fc£ for dis- 
order with fixed correlation parameter C, — cr/£ — 2, in 
d = 1,2,3 dimensions (top to bottom). For small or large 
universal limiting behavior is found, see text. The feature 
at ka — 1 is particular to the speckle disorder [Eq. (A6)]. 
Inset: the plot's path in the parameter space, Figure 4. 



The parameter space thus is a two-dimensional manifold. 
Nonetheless, it is useful to keep all three parameters to 
describe the various physical regimes. We have therefore 
found it convenient to map the entire parameter space 
to a hexagon, see Fig. 4. The three symmetry axes con- 
necting the vertices carry the three dimensionless para- 
meters, such that the extreme values and oo occur at 
opposite vertices. The symmetry axes perpendicular to 
the edges then represent the values of the length scales 
themselves, with their extreme values 0, oo taken on the 
entire edges. This construction is analogous to the rep- 
resentation of RGB color space by hue and saturation at 
fixed lightness. Each of the three dimensionless parame- 
ters describes one of the channels, e.g. fc£ the red channel 
with fc£ = mapped to cyan and fc£ = oo mapped to red. 
Similarly, ka and <r/£ define the green and blue channel, 
respectively, which completes the color space, as shown 
in Fig. 4. 



IV. RESULTS 

All of the physical quantities that we compute in the 
following depend crucially on the correlation parameter 
C = °7£ measuring the disorder potential correlation in 
units of the condensate healing length. In general, the 
results will even depend on the specific pair correlation 
function Cd(ka) defined in Eq. (71). For concreteness 
and direct applicability to cold-atom experiments, we will 
study in detail the case of optical speckle patterns, some 
properties of which are summarized in Sec. A. However, 
many analytical results in limiting cases are independent 
of the specific pair correlation and are thus universally 
applicable. 



A. Mean free paths 

1. Elastic scattering mean free path 

First, we evaluate the elastic scattering mean free path, 
in the dimensionless form l/(kl s ) = — 2lmS(k)/(hkv s ), 
which is the small parameter of weak-disorder expansions 
in standard quantum transport theories [59]. The only 
possibility for an imaginary part to occur in the on-shell 
self-energy (73) is by the imaginary part of the Green 
function, Im(e fe - e k+q + iOy 1 = -n5(e k - e k+q ), mul- 
tiplying the normal scattering amplitude. This restricts 
the integral over the intermediate state to the energy 
shell. There, the scattering element simplifies according 
to Eq. (36), and we find 

1 nv 2 k d a d f dn d t , li .» 2 „, nl . 9 , 

kf s = ^(TTIW J j^m,e?c d {2k*^). 

(81) 

Fig. 5 shows this inverse scattering mean free path 
plotted as function of fc£ for a speckle potential [Eq. (A6)] 
with a fixed correlation ratio C, = a/£ = 2 in dimensions 
d = 1,2,3. The fc-dependent fraction in front of the in- 
tegral in (81) ensures that the mean free path diverges 
both for very low and high momenta, where one recovers 
essentially a clean system. In-between, around k£ = 1, 
there is a minimum mean free path. Fig. 5 also shows a 
feature at ka = 1, which is specific to the speckle correla- 
tion function used here, due to the non-analyticity at the 
boundary of its support. Indeed, in d = 1, at this point 
the contribution of the backscattering process k — » —k 
becomes impossible at the level of the Born approxima- 
tion, which explains the kink at ka = 1. In higher di- 
mensions, the angular integration lifts the singularity to 
higher derivatives, so that it becomes less conspicuous. 
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Obviously, the mean free path (81) depends on the two 
dimensionless parameters fc£ and ka. Let us discuss some 
interesting limiting cases, which are easily identified in 
our parameter space representation. 

(i) First, the hydrodynamic limit £ — > is found on the 
lower right-hand edge of Fig. 4. Here, we recover exactly 
the elastic mean free path for scattering of sound waves 
[36], which only depends on the remaining parameter ka. 

(ii) Secondly, the limit £ — ¥ oo is found on the upper 
left-hand edge of Fig. 4. There, £ drops out together 
with the interaction energy gn from (81), which reduces 
exactly to the elastic mean free path for single-particle 
scattering, as calculated in [61, 62]. This reduced mean 
free path kl s can again only depend on ka. 

More generally, let us consider a fixed correlation ra- 
tio <t/£ of order unity, and look at the asymptotics as 
function of fc£, as plotted in Fig. 5. 

(iii) As fc£ — > 0, also ka — > 0. Then, the potential 
correlator becomes isotropic:, Cd(2fc<7 8in |) — > C t j(0), and 



pulls out of the integral over A(0, 9) 



integral yields J dfi^ cos 2 6 = Sj/d (a result showing the 
"concentration of measure" of the hypersphere's surface 
around its equator [63]). Thus, we obtain 



1 



2d(2n) 



j(k0 d ll + O(ktka)] 



(82) 



where S d is the unit sphere's surface (Si = 2, S% = 2n, 
S3 — 47r), and we define the effective (^-correlation disor- 
der strength 



v 2 s 



C d (0) 



a d V 2 
^(gnf 



(83) 



The low-fc behavior l/kl s oc (fc£) d indeed appears clearly 
in Fig. 5. In our parameter space Fig. 4, this is the 
asymptotic behavior of curves starting from the lower 
edge k = 0, for intermediate values of cr/£, i.e. rather in 
the center of the edge. Note that l^ 1 oc k^ 1 is propor- 
tional to the surface of the energy shell, i.e. the number of 
states available for elastic scattering. In the limit k — > 0, 
the elastic energy shell shrinks and the scattering mean 
free path diverges, even when measured in units of fc -1 . 

(iv) Conversely, as fc£ — > 00, also ka — > 00. But as 
soon as ka ^> 1, the disorder potential allows practically 
only forward scattering, and we can make a small-angle 
approximation to all functions of (9, such as 2ka sin | — > 
kaO. Then, the final result can be cast into the form 
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v 2 fdWO 

El (kaf 



[l + 0(l/k^l/ka)\ (84) 



where E a = h 2 /(ma 2 ) is the characteristic correlation 
energy [61, 62], and /i(C) = Ci(0) as well as 



[°° duu d - 2 (2( 2 -u 2 ) 2 



(27r) d -! (2( 2 + u 2 y 



(85) 



in d = 2, 3. And indeed, all curves in Fig. 5 show this 
decrease as fc -3 for large momenta. In our parameter 



space Fig. 4, this is the asymptotic behavior of curves 
arriving at the upper edge k = 00 for intermediate values 
of a/t- 

For extreme values of cx/£, i.e. towards the far left or 
right of the parameter space, one can also find interesting 
asymptotics when the wave length 2ir/k lies between a 
and £. 

(v) Consider first the case of a potential in the deep 
Thomas-Fermi regime and excitations with ka 3> 

1 3> /c£. This describes hydrodynamic excitations, i.e. a 
set of parameters approaching the right outermost vertex 
of Fig. 4 along the edge £ = 0. Then, by the same reason- 
ing as in the previous case (iv), strongly peaked forward 
scattering leads to l/kl s oc ka. This linear increase of the 
inverse scattering mean free path would naively predict 
infinitely strong scattering as k increases. Within the full 
description, however, this unphysical behavior stops as 
soon as kt; ss 1 is reached, crossing over to case (iv) with 
an even simpler description since now £ 3> 1 for which 



'0. This f d (oo) = Sa-i^ir) 1 -* duu d - 2 C d (u) hid = 2,3 



(vi) Conversely, consider finally a potential in the 
deep smoothing regime cr/£ -C 1 and excitations with 
ka <C 1 <C This describes particle excitations, i.e. 
a set of parameters approaching the left outermost ver- 
tex of Fig. 4 along the edge £ = 00. Then, by the same 
reasoning as in case (iii), isotropic scattering leads to 
l/kl B oc (ka) d ~ 4 . This low- A; divergence of the scattering 
mean free path naively predicts infinitely strong scatter- 
ing as ka <C 1 and severely limits the validity of simple 
perturbation theory for the single-particle case [61, 62]. 
Not so here, where the divergence is avoided once fc£ ~ 1 
is reached, and the interaction energy comes into play, 
crossing over to case (iii). 

In summary, our perturbation theory provides valid 
expressions for the elastic scattering rate or inverse mean 
free path in the full space of parameters. At a given 
value of C — the scattering rate is always a bounded 
function of fc, multiplied with the small parameter v 2 <C 
1, which vindicates the use of the momentum basis as a 
starting point for the perturbation theory. 



2. Transport mean free path 

If a disordered BEC is brought out of equilibrium, it 
will respond via its excitations. Therefore, it is of interest 
to study the transport properties of Bogoliubov excita- 
tions. In principle, a full-fledged quantum transport the- 
ory requires to calculate particle-hole propagators, which 
is certainly doable using the Hamiltonian (58), but be- 
yond the scope of the present article. Still, the previous 
results on the scattering mean free path can be gener- 
alized, with very limited additional effort, to the Boltz- 
mann transport mean free path [61, 62] that measures the 
diffusive randomization of the direction of motion. This 
transport mean free path l tr is defined by the same in- 
tegral expression (81), where the integrand is multiplied 
by a factor (1 — cost9). Fig. 6 shows a plot of I /kin as 
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Figure 6. (Color online) Inverse transport mean free path, 
(81) including the vertex correction factor (1 — cos9) under 
the integral, for a speckle disorder (A6) with fixed correlation 
parameter ( — er/£ — 2. 



function of kt; for a speckle disorder (A6) with fixed cor- 
relation parameter £ = cr/£ = 2 in dimensions d = 1, 2, 3. 

In one dimension, the only contribution to the inverse 
transport mean free path is the backscattering contribu- 
tion k — > — k, such that 



1 _ v 2 kad(2ka) 

Wr ~ T(l + fc 2 C 2 ) 2 



(86) 



Due to the finite support of the speckle correlation func- 
tion (A6a) backscattering is impossible for kcr > 1, and 
the inverse transport mean free path vanishes (within the 
Born approximation, and here we do not consider higher- 
order corrections to lt T [64, 65]), as clearly apparent from 
Fig. 6. In dimensions d > 2, there are finite contribu- 
tions from small scattering angles. Adapting the reason- 
ing of case (iv) from the previous section, one finds that 
1/fcZtr cx (fcc)~ 5 , with a prefactor that can be determined 
similarly. 



3. Localization length 

Just as phonons and particles, Bogoliubov excita- 
tions are expected to localize in disordered environments. 
Again, a full calculation is out of reach within the present 
article, but we can estimate the localization lengths of our 
Bogoliubov excitations in correlated disorder, based on 
general results on localization of particles and phonons. 

In one-dimensional disordered systems, the localization 
length l\ oc — 2Z tr , which describes exponential localiza- 
tion, is directly proportional to the backscattering length 
that we just calculated [ ]. From (86) we deduce 



1 



v 2 kcr d(2ka) 
kl^ c ~ T(l + fc 2 £ 2 ) 2 



(87) 



which agrees perfectly with [18], and also with [17], in 
the limits a — > and £ — > investigated there. Those 
phase-formalism approaches are particularly suited for 
ID systems, whereas our Green-function theory permits 
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Figure 7. (Color online) (a) Condensate density n(x)/n of a 
BEC (system size L m 1700£, periodic boundary conditions) 
in a blue-detuned speckle potential with amplitude V/gn = 
0.3 and correlation length a = £. (b)-(d) Selected Bogoliubov 
modes u v (x) (solid red) and v v (x) (dashed blue), obtained by 
exact diagonalization of the Bogoliubov-de Gennes equation 
(46). Low- and high-energy modes [(b) and (d)] are extended, 
while localization is most pronounced at intermediate energies 
[(c)] . The dashed gray line in panel (c) shows the exponential 
envelope predicted by Eq. (87). 



going to higher dimensions without conceptual difficul- 
ties. 

In two dimensions, the localization length is related to 
the transport mean free path via Zi oc = ^tr exp{^AZ t r}- 
This result can be derived using scaling theory argu- 
ments that hold very generally for single-particle exci- 
tations, and also the localization length of phonons has 
been shown to scale exponentially with / t r [67]. 

In three dimensions, localized and delocalized states 
can coexist, as function of energy separated by a mobil- 
ity edge. Phonons are localized at high energies and par- 
ticles are localized at low energies [67]. These opposite 
characteristics imply that when the disorder is increased, 
localized modes will start to appear at energies close to 
the point where Z tr is minimum. 

So in all dimensions, localization will be observed most 
readily, within finite systems, for modes that have the 
shortest localization length. Our results on the trans- 
port mean-free path in correlated potentials show (in 
agreement with the ID results of [18]) that modes around 
fc£ = 1 will be the first to appear localized. 

In one dimension, we have quantitatively verified the 
prediction (87) by means of an exact diagonalization of 
the inhomogeneous Bogoliubov-de Gennes equation (46), 
after solving the stationary GP equation (5) for the con- 
densate. Fig. 7 shows that indeed only Bogoliubov modes 
at intermediate energies fvjj v rs O.Ggn appear localized in 
the finite system. The observed localization length is 
compatible with the prediction, Eq. (87). 

All lengths calculated so far have the property that 
they diverge in the limit fc£ <C 1. In other words, sound 
waves can propagate over long distances and for long 
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times in these disordered systems. It is therefore mean- 
ingful to compute the renormalized speed of sound. 



B. Speed of sound 

As shown in Sec. HID, the disorder potential shifts 
the dispersion relation by Aejt = ReS(fc). Using Eqs. 
(77) and (73), the relative dispersion shift takes the form 



Ae fc 



d d q 
(2ir) d 



ZkqC d (qa). 



The kernel z kq obtains from the real part of the on-shell 
kernel Z k{k+q)k (t k ), Eq. (72): 
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Cfe — e k+q (-k + (k+q 



, (2) 

' J k(k+q)k 



(89) 



with the envelopes defined in Eqs. (35a), (35b), and (42). 
P denotes the principal value. These equations reduce to 
much simpler expressions in different limiting regions of 
the parameter space, Fig. 4. We mostly focus on the low- 
energy, sound excitations that are of primary interest. 
Analytical results will be confronted with data from a 
numerical simulation in Sec. IV B 2 below. At last, we 
show in Sec. IV B 4 that our theory also covers the case 
of weak lattice potentials. 



1. Limiting cases 

Similar to the proceeding in Sec. IV A, we compute an- 
alytical results in the limiting cases located at the edges 
and corners of the parameter space, Fig. 4. 

(i) We start with the hydrodynamic limit £ — > 0, the 
lower right edge of Fig. 4. The kernel (89) simplifies to 



1 (k 2 + k-q) 2 
2k 2 q 2 + 2k-q ' 



(90) 



after which Eq. (88) reproduces exactly Eq. (29) of Ref. 
[36]. Notably, the dispersion shift is negative in all di- 
mensions and for any value of cr/£ 3> 1, as anticipated in 
the schematic plot of Fig. 2(a). The limiting values are 



Ae fc 



-l/(2d), 
-(2 + d)/8, 



fc£ < ka < 1 
fc£ < K ka 



(91a) 
(91b) 



These limiting values are expected to hold over an ex- 
tended range of ka. Thus, they define a shift in the 
local slope of the dispersion relation. In other words, 
fc-modes in that particular range have a modified sound 
velocity. The magnitude of the correction depends sig- 
nificantly on the excitation's ability to resolve the corre- 
lations (ka 3> 1) or not (ka < 1). As noted in [36], in 
the latter case the correction decreases with dimension, 



but increases in the former, implying that the curves for 
different dimensions must cross around ka = 1. 

In passing, we stress that even in the very long-range 
correlated limit <r/£ — > oo, these results are not trivial. 
Indeed, one could try and use a simple static local density 
approximation (LDA) in order to derive the result (91b) 
for correlation lengths much longer than the excitation 
wave length. In this regime, the background appears 
locally homogeneo us to the wave, and the local sound 
velocity c(r) = \J gn(r) jra is proportional to the con- 
densate field amplitude $(r). Thus, LDA expects Ac/c 
to be given by < J > / < I > , which can be easily computed from 
Eq. (10) to yield Ac L da/c = -v 2 /8. But this fails to 
reproduce Eq. (91b). Indeed, static LDA cannot capture 
the scattering dynamics (shown by the first diagram of 
Eq. (69)), which is essential for correctly determining the 
sound velocity. 

(ii) In the regime of particle-like excitations fc£ — > oo 
(covering the cases (ii), (iv), and (vi) of Sec. IV A), the 
Hamiltonian (1) becomes non-interacting. Consequently, 
one should expect the entire Bogoliubov problem to re- 
duce to the problem of single particles in disorder. In- 
deed, Bogoliubov excitations in the particle regime see 
both, the external potential and the condensate back- 
ground. Sampled at high wave numbers k£ 3> 1, the con- 
densate background is smooth and cannot induce scat- 
tering. Fittingly, we found in Sec. IV A 1 that the elastic 
scattering mean free path reduces in this limit to the 
single-particle expression. In contrast, the deeply inelas- 
tic processes contributing to Eq. (88) remain sensitive to 
the condensate background, as encoded by the anoma- 
lous and second-order couplings, (35b) and (42), which 
do not simply vanish in the limit £ — > oo. We find that 
the leading-order correction to the dispersion relation (for 
fc£ S> 1 and ka not too small), 
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(92) 



is independent of k. Incidentally, it is exactly opposite to 
the negative average shift //( 2 ) of the chemical potential, 
from Eq. (15), for fixed average density. At fixed chemical 
potential, the dispersion shift (92) would even be twice as 
big. Note that this shift cannot be naively accounted for 
by an overall shift gn — > gn + in the clean dispersion, 
Eq. (25). Just like the wrong LDA attempt to explain the 
sound velocity, discussed above, such a reasoning misses 
the essential scattering dynamics. 

Eq. (92) differs also from the chemical potential shift 
for noninteracting particles in disorder [55]. But it must 
be kept in mind that the disorder expansion of the Bo- 
goliubov Hamiltonian in Sec. II B 3 was performed under 
the assumption V <C gn. For single particles, the inter- 
action energy gn goes to zero, i.e. the ratio of V and gn 
would have to be reversed. Therefore, our perturbative 
theory cannot be expected to apply universally in this 
regime. 

In any case, the main effect of disorder in the single- 
particle regime is to yield the finite scattering rate cal- 
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Figure 8. (Color online) Relative correction (93) of sound 
velocity due to a speckle disorder potential (A6) as function 
of the correlation ratio £ = cr/£ in d = 1,2,3. The exact 
formulae are given in Eq. (All). Dashed: universal limits, 
collected in Table I. Inset: same data around the origin, 
showing the rapid departure from the leading-order estimate 
[14, 19, 33]. 



culated in Sec. IV A, but it only produces a very small 
shift in dispersion. In addition, these high-energy excita- 
tions are less important for low-temperature properties 
of BECs, and will not be considered in the remainder of 
this work. 

(iii) Let us turn to the sound-wave regime fc£ <C 1. In 
case (i) and Ref. [36], this has been achieved by sending 
the healing length £ to zero, thus yielding the dispersion 
as function of ka, but only for rather long-range corre- 
lated potentials with a 3> £. In order to cover arbitrary 
correlation ratios £ = we now change the point of 
view and take k — > 0. This allows in particular to reach 
the case a <C £,fc -1 of truly ^-correlated disorder that 
was inaccessible to [36]. The kernel (89) simplifies, and 
we find the relative shift in the speed of sound 
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C d (qa). 
(93) 



P = Z(fe, q) is the angle between the direction of propa- 
gation and q. In contrast to the hydrodynamic case (i), 
there are now two competing contributions with opposite 
sign (for an interpretation of these contributions, cf. the 
end of Sec. IV B 4 below). In case of isotropic correlation, 
the angular integral maps cos 2 f) to 1/d. Then, only for 
d = 1 is the radial integrand strictly negative, and Ac is 
negative as well. For d > 1, the radial integrand has no 
definite sign. 

To survey the possible outcomes, we plot in Fig. 8 the 
correction (93) to the speed of sound caused by isotropic 
speckle disorder, Eq. (A6), as function of the correlation 
ratio £ = c/£. The curves can actually be given in closed 
form, see (All), but the details depend of course on the 
specific correlation. In contrast, we can extract universal 
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Table I. Limiting corrections of the speed of sound, corre- 
sponding to the dashed limits of Fig. 8. For uncorrelated 
disorder, the correction is proportional to v 2 of Eq. (83). 



limits for very small or very large £. 

In the long-range correlated limit £ — ► oo, found on 
the right edge of the plot, the correlator Cd acts as a 
^-distribution, which leads to Ac/c = — v 2 /(2d). This 
value coincides with the hydrodynamic limit (91a), as it 
should. 

In the opposite limit ( < 1 of <5-correlated disorder, 
the correlator Cd{qcr) — > C<z(0) can be pulled out of the 
integral, which contributes a numerical prefactor to the 
expected scaling with the disorder strength v 2 defined 
in Eq. (83). These results are plotted as dashed lines 
in Fig. 8 and collected in Table I. Again, these results 
cannot be found by LDA, which would have to assume 
that the system is homogeneous on relevant length scales 
(o 3> £, 2n/k), an assumption that is always violated by 
the sound-wave limit 2?t /k oo. 

Our result for a <C £ in d = 3 reproduces the value 
known from Refs. [14, 19, 33]. Interestingly, this is the 
only case where the correction to the speed of sound is 
positive. Actually, this particular numerical value is of 
rather limited use since the Taylor expansion at the ori- 
gin is converging very slowly, and already a minor cor- 
relation can make a major difference, as shown by the 
inset in Fig. 8. Our results, Eqs. (88) and (93), hold for 
a much larger range of parameter values and arbitrary 
dimensions, which accomplishes one of the main goals of 
this work. 



2. Numerical mean-field study of the sound velocity 



We confront the theoretical predictions (88) and (93) 
with data obtained by a numerical simulation in d = 1 
on the mean-field level, using the time-dependent Gross- 
Pitaevskii equation. This numerical calculation consti- 
tutes an independent check since it relies neither on the 
linearization in the excitations, nor on perturbation the- 
ory in the disorder potential, which are the two approxi- 
mations of our analytical theory. 

The numerical procedure has been briefly described in 
Ref. [36]. We generate a ID speckle disorder potential 
with correlation length a by Fourier transformation from 
the set of random complex field amplitudes [see Eq. (A5)] 
with all k < cr _1 . Then the condensate ground state <&(x) 
solving the GP equation (5) is computed by imaginary- 
time evolution, using the fourth-order Runge-Kutta algo- 
rithm, while keeping the wave function normalized [46]. 
Onto this disordered ground state, a plane-wave Bogoli- 
ubov excitation is superposed, with a small, but finite k 
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Figure 9. (Color online) Relative correction of the Bogoli- 
ubov excitation dispersion relation due to ID speckle disorder. 
The full formula (88) for fc£ = 0.05 [solid black] crosses over 
from the limiting case (iii) of low-energy excitations [dashed 
green, Eq. (Alia)] to the limiting case (i) of the hydrody- 
namic regime [dotted blue, Eq. (A7)]. The inset shows the 
corresponding trajectory in parameter space. The numerical 
results (cf. Sec IV B 2) for blue- and red-detuned speckle with 
v = +0.03 [blue lines] and v = —0.03 [red triangles], agree 
fully with the analytical theory. 

and amplitude T. In cold-atom experiments, such an ex- 
citation is routinely imprinted using Bragg spectroscopy 
[68-71] . The Bogoliubov transformation (20) requires the 
imprints in density and phase to be 

r 

Sn(x) = 2\/na,k rcos(fca;), S(p(x) = —= — sin(fcr). 

(94) 

In the sound-wave regime where si < 1 [cf. Eq. (23)], the 
phase modulation has a much larger amplitude than the 
density modulation, and we choose T = 0.3y/nakV. Then, 
the real-time evolution under the GP equation is com- 
puted using again the fourth-order Runge-Kutta algo- 
rithm. The excitation propagates, with a modified speed 
of sound, surviving over a long course of time given by 
the inverse elastic scattering rate 7^7 1 [cf. Eq. (76)]. In 
order to extract the eigenenergy e/., the deviation <5$(x, t) 
from the ground state is translated into Bogoliubov ex- 
citations by means of Eq. (20). Monitoring the phase of 
7^ cx e - ie kt/h over time, we extract the phase velocity 
Vph = £k/frk by linear regression. Then, the procedure 
is repeated for different realizations of disorder with the 
same correlation length tr, leading to a whole distribution 
of values, from which we compute the average Ac/c. As 
shown in Ref. [36], for weak disorder the distribution is 
clearly single-peaked and allows for meaningful averages. 

Figure 9 shows the numerical data on top of the the- 
oretical predictions, as function of a. Since fc£ = 0.05 is 
fixed, the curve can be read as a function of ka at very 
small £ (near the hydrodynamic limit (i) of above) or as a 
function of cr/£ at very small k (near the low-energy limit 
(iii) of above). The inset of Fig. 9 depicts the correspond- 
ing trajectory in parameter space. The full prediction 
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Figure 10. (Color online) Relative correction (88) of the dis- 
persion relation Aek for different correlation ratios £ = <r/£ = 
0.2, 0.5, 1, 2, 5, 10, 20 (top to bottom), for ID speckle disorder, 
Eq. (A6a). The results of Figure 8 appear at the edge fc£ = 
(circles), whereas the results shown in Figure 9 are found at 
fc£ = 0.05 (triangles). Around the points ka = 1 (diamonds), 
the correction behaves non-monotonically. The points with 
errorbars close to the curve £ = 2 show data from the exact 
diagonalization of the Bogoliubov-de Gennes eq. (46) (system 
size L « 300£, correlation a — 2£) for blue- as well as red- 
detuned speckle disorder. Each point represents the energy 
shift of the two modes v — 2j — 1, 2j corresponding to the 
degenerate modes kj = ±2irj/L of the homogeneous system. 
The data has been averaged over a large number r of realiza- 
tions (rL/a ~ 1.9 x 10 4 ). Errorbars show the estimated error 
of the mean value. 

(88) for fc£ = 0.05 is plotted as a black line. The limiting 
cases are available in closed form: Eq. (A7) for £ = as 
function of ka and Eq. (Alia) for k = as function of 
cr/£. Independently of the details, we find of course the 
relevant universal limits of Sec. IV B 1 above. 

The numerical data, shown for v = +0.03 (blue 
straight marks) and v = —0.03 (red triangular marks), 
follow the analytical prediction very well. Interestingly, 
the data points of red and blue detuning, with opposite 
sign of v, tend to lie on opposite sides of the curve, in- 
dicating beyond-Born effects of odd order v 3 , which are 
expected for a speckle potential with its asymmetric on- 
site distribution (A2). Attentive readers will also notice 
that the data points are shifted asymmetrically with re- 
spect to the curve, which is an effect of order v 4 . 

3. Disorder-shift of Bogoliubov spectrum 

Finally, we explore the correction of the dispersion re- 
lation (88) as function of k£. Figure 10 shows a family of 
curves for different correlation parameters £ = <t/£ in one 
dimension. This plot actually contains the information 
of the previous Figs. 8 and 9, which are taken at fixed 
values of fc£ = and fc£ = 0.05, respectively. 

The disorder correction passes through a non- 
monotonic feature around ka — 1 (marked by open dia- 
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monds), and finally diminishes with increasing 

In dimensions d > 1, the curves have a slightly differ- 
ent shape. The curve with £ = 20, for example, starts 
with the limit Eq. (91a) at fc£ — > and passes through 
Eq. (91b) at fc£ « 0.15. Thus, while the ID curve starts 
at —1/2 and passes through —3/8, the corresponding 
curve in d = 3 will start at —1/6 and pass through an 
extremum around —5/8, before diminishing in the parti- 
cle regime. Also, the sharp speckle features get washed 
out in higher dimensions. 

Complementary to the previous numerical study at 
constant kt;, we can verify our predictions also by an 
exact diagonalization of the one-dimensional disordered 
Bogoliubov-de Gennes equation (46). Thus, we obtain 
the whole spectrum of the system characterized by its 
correlation ratio £ = cr/t;. For the case C = 2, Fig. 10 
shows excellent agreement between prediction and nu- 
merics. As observed previously, the data for red- and 
blue-detuned speckle lie to both sides of the 0(v' 2 ) pre- 
diction, indicating effects of odd orders. 

The excellent agreement between both numerical ap- 
proaches and analytics demonstrates that our perturba- 
tion theory to order v 2 gives an impressive account of all 
relevant effects over the entire parameter space we set 
out to cover. 



for 

V 2 a d C d (qa) = { ^f- ]T V 2 [S(q - K 3 ) + S(g + K s )] . 

3=1 

In the sound-wave limit k 0, the speed-of-sound cor- 
rection (93) thus reads 

Ac 'j K 2 e cos 2 f3j ] V 2 

c j^\{2 + K 2 ef (2 + K 2 e) 2 j (gn) 2 1 ' 

where fij is the angle between the propagation direction 
and the lattice direction e.j. In the case V2 = V3 = 0, 
the sound velocity along the x\ direction reproduces the 
result of Ref. [37, Eq. (52)]. Also, Eq. (97) is precisely 
Eq. (22) in Ref. [38] for the potential (96). This work 
of Liang et al. [38], while equivalent to ours as far as 
the perturbative approach is concerned, has the nice fea- 
ture that it permits to attach a thermodynamic meaning 
to the two antagonistic contributions appearing in Eq. 
(97): the first, positive term stems from the disorder- 
shift of the compressibility re, while the second, negative 
one goes back to the change in the effective mass m*. 
Together, these quantities determine the speed of sound 
c = l/V rem*. 1 



4- Weak lattice potentials 



C. Average density of states 



The inhomogeneous Bogoliubov Hamiltonian (58) ap- 
plies to arbitrary external potentials. In particular, it 
covers the important class of weak optical lattice poten- 
tials that was studied in [37, 38]. In a sense, understand- 
ing the lattice is a first step to understanding disorder, 
which can be seen, by virtue of Fourier decomposition, as 
a superposition of lattices with suitably chosen random 
amplitudes and phases. 

To substantiate this connection, we briefly show that 
our formalism reproduces the results of [37, 38] for the 
speed of sound in weak lattices. In d dimensions, a sep- 
arable lattice potential V(r) = Ylj=\ Vj C0S (Kj x j) with 
wave vectors Kj = Kj&j has Fourier components 



1 d 



(95) 



The whole formalism developed in Sees. II and III applies 
also to this potential. Notably, the equation of motion 
(63) is still solved by the perturbative expansion (67), 
now without the need for averaging over disorder. The 
periodic potential scatters Bogoliubov excitations elasti- 
cally only at the edges of the Brillouin zone, kj = ±Kj. 
Away from the edges, the dispersion relation is again de- 
termined by the diagonal elements of the Green function 
G(e)fcfc. Expressions like (73) and (88) are still valid, 
where the correlator should now be understood to stand 



Lastly, we turn to the average density of states (AV- 
DOS), Eq. (78), at very low energies or k — > 0. As shown 
above, the dispersion is linear in this limit, so that the 
AVDOS p(e) = J ^S(e-hck) = p{e){c/c) d necessarily 
has the lowest-order correction 



Ap(0) = Ac 
P (0) 



(98) 



In other words, a reduced sound velocity entails an en- 
hanced density of states and vice versa, which is the obvi- 
ous conclusion one can already draw from the schematic 
plot anticipated in Fig. 2. It is instructive to check that 
identity (98) also follows from the general equation (79): 
There, the linear dispersion implies kd^k = £fc, and the 
fraction approaches Alk / e& . When acting on this regular 
function, the operator kdk inside the bracket evaluates to 
zero at k = 0, and we arrive at (98). Thus, the AVDOS 
shift is entirely determined by the sound velocity shift 
(93); see also the analytic solutions (All). 

In the hydrodynamic regime fc£ <C 1 realized by £ — s- 
at finite ka, again kdk^k — £fe = hck in Eq. (79), which 
then reproduces our previous results [36]. The limiting 



1 The compressibility correction in Ref. [38] appears with a wrong 
sign in Eq. (CI), but (C6) agrees with our result. 
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values for small or large energy compared to the hydro- 
dynamic correlation energy he/ a, corresponding to (91), 
are 



Ap(e) 



1(2 



d), 



e <C he/ a, 
e 3> he/ a. 



(99) 



In between these limits, the correction Ap/ p as function 
of e can show a surprisingly rich behavior, depending on 
the potential correlations. In three dimensions, the cor- 
rection is smooth and monotonic, but in one dimension, 
speckle correlations are responsible for an unexpected, 
non- monotonic behavior with a sharp feature at k e a = 1, 
as discussed in [36]. 

Note that these simple perturbative results are indeed 
expected to hold at low energy or fc£ <C 1, in stark con- 
trast to the case of single particles in disorder, where the 
DOS has a (nonperturbative) Lifshitz tail at low energy 
[72, 73]. In the present case, the interparticle repulsion 
screens the disorder very effectively at low energy, such 
that localization effects are absent (cf. the diverging lo- 
calization length of Sec. IV A 3) [12, 73]. A transition 
to the Bose glass phase [1, 2] occurs only for stronger 
disorder or much weaker interaction, where the Bogoli- 
ubov theory developed here breaks down, and different 
approaches are needed [74-76]. 



V. CONCLUSIONS AND OUTLOOK 

In conclusion, we have formulated a comprehensive Bo- 
goliubov theory of inhomogeneous Bose-Einstein conden- 
sates. This analytical theory describes the elementary 
excitations of condensates with s-wave interaction, de- 
formed by weak external potentials with arbitrary spatial 
correlations and in arbitrary spatial dimension. Expand- 
ing the many-body Hamiltonian around the deformed 
ground state, we have obtained the inhomogeneous Bo- 
goliubov Hamiltonian. We have justified our choice of 
the basis of density and phase fluctuations that ensures 
proper orthogonality between the excitations and the in- 
homogeneous Bogoliubov vacuum. Expressed in terms 
of Nambu-Bogoliubov spinors, all effects of the external 
potential can be collected into a scattering vertex that is 
non-perturbative in the external potential. A fully ana- 
lytical formulation has been achieved up to second order 
in weak potentials, allowing in principle an expansion to 
even higher orders. 

From this fundamental Hamiltonian, one can derive 
numerous physically relevant quantities by means of stan- 
dard perturbation theory. This paper has been devoted 
to a detailed discussion of the single-excitation disper- 
sion relation. We have calculated the mean-free path and 
renormalized speed of sound together with the resulting 
average density of states, over the full parameter space of 
the disordered Bogoliubov problem, with numerous ana- 
lytical results in limiting cases. It turns out that the fre- 
quently investigated case of (5-correlated disorder in three 



dimensions, with its positive shift in the speed of sound, 
is far from generic. Over most of the parameter space, 
the speed of sound is reduced. We have confirmed these 
predictions in detail by mean-field numerical simulations 
as well as exact diagonalization for the experimentally 
relevant case of correlated speckle disorder in d = 1. 

Strictly speaking, the present work is incomplete with- 
out proving that the weak disorder under consideration 
causes only a small condensate depletion. Indeed, the Bo- 
goliubov ansatz (3) relies on the macroscopic population 
of the condensate mode. In a pure mean-field description, 
all atoms are in the condensate at zero temperature. But 
due to the effect of interaction, even at zero temperature, 
there is a finite fraction of particles not in the conden- 
sate, which constitutes the so-called quantum depletion 
that can be calculated within Bogoliubov theory [31, 32]. 
The quantum depletion should be small, thus providing 
an important, self-consistent check of the theory's valid- 
ity. 

In the homogeneous setting, the mean-field conden- 
sate forms in the homogeneous mode, i.e., the zero- 
momentum state. Thus, the quantum depletion consists 
of all particles with non-zero momentum. The density of 
these particles can be easily calculated within Bogoliubov 
theory [31, 32], with the result (we take d = 3 here) that 
the fractional quantum depletion 8n/n — 8(nag) 1 / 2 /3- v /7r 
is proportional to the root of the gas parameter and thus 
very small, especially so for dilute and weakly interacting 
cold gases. 

It is not immediately obvious how to generalize the 
recipe "count all particles with non-zero momentum" to 
the inhomogeneous case. The vast majority of works ded- 
icated to the inhomogeneous Bogoliubov problem simply 
calculates the same quantity, namely the number of par- 
ticles with non-zero momentum. But one has no means of 
knowing whether these particles belong to the deformed 
condensate or to the true, disorder-induced quantum de- 
pletion. And really, the supposed "depletion" calculated 
by Huang and Meng, followed by [13-16], involves only 
the condensate deformation at fixed chemical potential, 
which is a mean-field effect as described in Sec. II A. Only 
few authors seem to have clearly recognized that this sup- 
posed depletion describes merely the non-uniform density 
of the condensate [33] . 

In order to assess the true condensate depletion, one 
has to determine the density of particles not in the con- 
densate at all, irrespective of their particular momentum. 
With the general Bogoliubov Hamiltonian at hand, we 
have calculated this disorder-induced quantum depletion 
[ ] . We find that it is much smaller than the mean-field 
condensate deformation. This is no surprise: To first 
order, the external potential merely deforms the conden- 
sate. The scattering of particles out of the condensate 
is a second-order effect, mediated by the interaction be- 
tween particles and the condensate. Details of the full 
calculation, including finite-temperature effects, will be 
discussed in a forthcoming publication [78]. 

These results validate and strengthen the Bogoliubov 
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approach, and we expect that the theory we have devel- 
oped here should fare very well in describing the excita- 
tions of inhomogeneous BECs. As an immediate exten- 
sion of the present work, finite-temperature effects can be 
captured very straightforwardly by the Matsubara for- 
malism [55], allowing the calculation of the heat capacity 
and many other (thermo-)dynamic response functions. 
This is left for future work. 
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Appendix A: Optical speckle disorder 
1. Statistical properties 

The optical speckle field of a laser defines a disorder 
potential with very well controlled statistical properties 
[24, 64, 79, 80]. When coherent laser light is directed on 
a rough surface or through a diffusor plate, the elemen- 
tary waves originating from different points have random 
phases and form a random interference pattern in the far 
field. By virtue of the central limit theorem, the result- 
ing field £(r) is a complex Gaussian random process. 
For notational simplicity, we consider a scalar field with 
dipole coupling to a single atomic transition and neglect 
polarization issues. The electronic atomic ground state 
is then subject to the light-shift potential induced by the 
intensity I(r) = |£(r)| 2 [81]. Magnitude and sign of this 
potential depend on the laser detuning from the dipole 
transition frequency. For a far-detuned potential, one has 



V(r) = V 



I(r) 



(Al) 



The prefactor V contains the atomic polarizability, be- 
sides all other proportionality factors, and we have 
shifted the potential to zero average, V(r) = 0. The 
magnitude of V gives the variance of the potential fluc- 
tuations, V 2 = V(r) 2 , that can be readily adjusted in 
the experiment by changing the overall laser intensity. 
Because the intensity is the modulus square of the Gaus- 
sian field £, the potential has a skewed on-site probability 



distribution for w — V(r)/V, 

P(w)dw = Q(w + l)e- (w+1) dw. 



(A2) 



A blue-detuned potential with V > consists of high 
potential bumps rising over a flat baseline. A red-detuned 
potential with V < rather corresponds to a random set 
of deep wells. 

The potential correlation between different spatial 
points is captured by the pair correlator 



C d (r/a) = V(r)V(0)/V 2 



(A3) 



that decays from the starting value Cd{0) = 1 to zero on 
the scale of the correlation length a. In momentum rep- 
resentation, the normalization implies the useful identity 



W d(<z<7) 



1. 



(A4) 



A speckle pattern's correlation is entirely determined 
by the fact that the Fourier components of the field are 
independent, complex Gaussian random variables with a 
pair correlation 



>y(k)6 



kk' 



(A5) 



that defines the degree of coherence j(k). In one dimen- 
sion, for a rectangular source, the degree of coherence 
j(k) — 70(1 — fcer) [62] has uniform weight inside the 
interval of allowed fc-values. The correlation length a de- 
pends on the laser wave length and the imaging optics 
and is typically of the order of 1 /mi or smaller [82] . It is 
not easy to create an isotropic multi-dimensional speckle 
field in the laboratory [26] . For the purpose of the present 
paper, we follow Pilati et al. [75] and define the isotropic 
speckle disorder on the grounds of Eq. (A5) with the same 
j(k) in all dimensions. This definition grasps the essen- 
tial feature of speckle disorder, namely the finite support 
of its power spectrum. The extension to a more realis- 
tic, possibly anisotropic correlation function for a given 
experimental configuration is straightforward within the 
present formalism. 

The potential correlator Cd(ka), used in Eq. (71) and 
the following, then obtains as the auto-convolution of the 
field correlator j(k), i.e. of a d-dimensional ball of radius 
a -1 . Thus, Cd{qu) is centered at q = and vanishes like 
(2 — qa)^~ at q = 2er~ 1 . Explicitly, in d = 1,2,3 one 
has 



C 1 (q<T) = -(2-qa)Q(2-qa), 



C 2 {qa) = 
C 3 (qa) = 



8arccos (^ 



2qay/A- q 2 a 2 



(A6a) 

6(2 -qa), 
(A6b) 



Sir 2 



[2-qaf{A + qa)®{2-qa). 



(A6c) 



Figure 11 shows these isotropic correlation functions, 
multiplied with the d-dimensional integration element 
Sd{q°') d ~ 1 1 '(2'?r) d , where Sd is the unit-sphere surface 
(Si =2,S 2 = 2tt, S 3 = 4tt). 
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(27T) 




Figure 11. Isotropic speckle correlation function as function 
of reduced momentum in d = 1, 2, 3. In order to show compa- 
rable scales, we plot [S d /(2n) d ](qa) d - 1 C d (qa), including the 
d-dimensional volume integration element. 



2. Analytical dispersion corrections 

The speckle correlation functions (A6) yield closed- 
form expressions for the dispersion correction (88) in im- 
portant limiting cases. 



Hydrodynamic limit £ = 



non-analyticity of (A8) at ka = 1. 

In higher dimensions, the integral (88) gets more com- 
plicated, but we find partial analytical results in two di- 
mensions. Denote the angular part of the integral (88) 
over the hydrodynamic kernel (90) by ^2(9)- We drop 
the principal value P in Eq. (90) and re-insert the in- 
finitesimal imaginary shift in the denominator. Then, we 
can compute the angular integral analytically as a closed- 
path integral in the complex plane z = e 1 ^ , (i — Z(fc, q) 



A 2 (q) 



1 So 



2 (2tt) 5 



f-V 

\2k) 



■ l~g 2 /(2fc 2 ) 
q V<Z 2 - Wf 



(A9) 

In general, the last term is too complicated for the re- 
maining radial integral to be evaluated in closed form. 
For ka > 1, however, the integrand of Eq. (88) is re- 
stricted to q < 2/ a < 2k. There, the last term in (A9) is 
imaginary and does not contribute to ReE, such that 



Ae fc 



v 



1 - 



1 



ka > 1. 



(A10) 



In the AVDOS (79), this leads to a totally flat plateau, 
as shown in Fig. 5 of Ref. [36] . 



In d = 1, the correlation function (A6a) is piecewise 
linear, and the integral (88) over the hydrodynamic ker- 
nel (90) yields the dispersion shift [36, 83] 



In 



1 — ka 



1 + ka 



2^2 



k 2 a 



In 



1 - k 2 a 2 



k 2 a 2 



(A7) 



This result is plotted as the blue dotted curve in Figure 9. 
It is non-analytic at ka = 1, corresponding to the non- 
analyticity of the speckle pair correlation function (A6a) 
at the boundary of its support. Using Eq. (A7), we can 
also write down the AVDOS shift (79) in closed form: 



Ap(ehc/a) 
p(ehc/a)v 2 



In 



1 -£ 
1 + e 



3^ 

8 



In 



1-e 2 

c 2 



(A8) 



It shows a pronounced dip around e ~ 0.7 and a sharp 
logarithmic divergence at e = 1 [36], resulting from the 



b. Lowest-energy excitations, k — > 

Also in the low-energy limit (iii) of Sec. IV B 1, we 
can find analytical solutions. The integral (93) with the 
speckle correlator (A6) evaluates to closed form in all 
relevant dimensions: 



Ac 
cv 2 
Ac 
cv 2 
Ac 
cv 2 



cot" 1 (z) 



1 z 
31 + z 2 
,2zy/l + z 2 - 1 - 2z 2 



VTTz 2 
5cot _1 (z) 
2z 



- (6 + 7z 2 )log 



d = 1 (Alia) 

d = 2 (Allb) 

1 + z 2 " 



d = 3 (Allc) 



with z = C/V2- All three cases are plotted in Figure 8; 
the result (Alia) also features in Figure 9. 
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